%% 人工神经网络Artificial neural network
%% 1.全连接神经网络Fully connected neural network
%%  1.1.MP模型简介
% 本节内容衔接智能预测与分类理论部分。
% 随机优化算法被广泛应用在机器学习，人工智能领域。本节用sign函数作为激活函数
% 模拟单个人工神经元建立MP模型，分析不同的随机优化算法。
clear,clc
% 设有原始数据：
n=1200;% 观测次数
m=4;% 自变量个数
x=-1+2*rand(n,m);% 观测数据为m项，每项n个数据
y=zeros(n,1);
w=[-1,-1,1,1];b=0;% 预设w,b设计原始数据
for i=1:n
    y0=sum(w.*x(i,:))+b;
    if y0<0
        y(i)=-1;% 对数据预分类
    else
        y(i)=1;
    end
end
x=[ones(n,1),x];
f=@(w,x)sign(sum(w.*x));% 简单MP模型，一般情况下可以进行二分类计算。
% 此处定义损失函数：
L=@(w,x,y)1/4*(y-f(w,x))^2;% 当分类结果正确时计算为0，否则计算为1
% 利用损失函数梯度下降设计学习算法:
gLw=@(w,x,y)1/2*(f(w,x)-y)*x;% 梯度
% 设置初始学习率(即步长)
alpha=0.001;
mn=1e3;% 总迭代次数
w=rand(1,m+1);w=w/norm(w);% 初始化w
accu=zeros(mn,1);% 统计每次迭代准确率
for i=1:mn
    tic
    Lpa=0;
    for i1=1:n
        y_pre=f(w,x(i1,:));
        Lp=L(w,x(i1,:),y(i1));
        dL=gLw(w,x(i1,:),y(i1));
        w=w-alpha*dL;
        w=w/norm(w);
        Lpa=Lpa+Lp;
    end
    accu(i)=(n-Lpa)/n;
    toc
end
% 由以上结果分析可知，梯度下降算法对于严格线性可分数据有较好的收敛速度和准确
% 率，但收敛算法容易受噪声影响。一般情况下，当数据有噪声时可以采用随机梯度下
% 降法，在样本中随机抽取一组数，计算该组数梯度平均值作为梯度下降方向进行迭代。
%梯度下降的拓展策略：
% 在掌握基本的梯度下降迭代方法后，可以采用一定策略对迭代进行加速或使迭代后参
% 数更加准确。
% (1)动量法：动量法迭代格式为：v=miu*v-alpha*dL;x=x+v;动量法引入了速度变量v
% 它使参数w的移动具有的一定的方向和速度，当迭代方向基本不变时，参数会加速变
% 化，从而减小迭代次数。
% (2)AdaGrad：迭代格式为：x=x-alpha./(sqrt(G)+eps)*dL;G=G+dL.^2;可以看到，
% AdaGrad方法对每一个参数迭代学习率(步长)进行了变化，当该参数对应变量偏导过
% 大时适当减小该参数对应学习率；反之则反。
% (3)RMSProp：迭代格式为：x=x-alpha./(sqrt(G)+eps)*dL;G=rho*G+(1-rho)*dL.^2;
% RMSProp是对AdaGrad适当的修改，防止出现训练后期步长过小的情况。
% (4)AdaDelta:迭代格式为：G=rho*G+(1-rho)*dL.^2;deltax=-(sqrt(D)+eps)/(sqrt(G)+eps)*dL;
% D=rho*D+(1-rho)*deltax.^2;x=x+deltax;AdaDelta算法改善了AdaGrad算法步长单调
% 下降的缺陷。
% (5)Adam：Adam算法本质上是带动量项的RMSProp算法，它利用梯度的一阶矩和二阶矩
% 估计动态调整每个参数的步长。令初值s=0,G=0,有迭代格式：
% s=rho1*s+(1-rho1)*dL;G=rho2*G+(1-rho2)*dL.^2;s=s/(1-rho1);G=G/(1-rho2);
% x=x-alpha./(sqrt(G)+eps)*s;
%%  1.2.全连接神经网络
% 本节利用全连接神经网络进行回归分析。全连接神经网络基本单元是MP模型，由3层组
% 成，分别为输入层，隐含层，输出层，其中隐含层可以为多层，隐含层和输出层具有
% 神经元。
% 下方通过数据学习训练过程说明全连接神经网络运算过程。
clear,clc
% 设有观测数据：
n=1500;% 观测次数
m=10;% 自变量个数
x=rand(n,m);% 观测数据为m项，每项n个数据
y=zeros(n,m);
for i=1:m
    y0=i^2*x(:,i)+2*i;
    y(:,i)=y0/norm(y0);
end
y=sum(y,2);
y=1./(1+exp(-y));% 规范化y
x=[ones(n,1),x];% 添加偏置项
% 设计神经网络模型：
num_net=struct('dim',m, ...% 创建结构体，声明输入向量维度为m
    'num_hid1',1,'num_hid2',8,'num_out',1);% 隐含层层数为1，隐含层神经元个
% 数为2，输出个数为1.
w1=rand(num_net.num_hid2,m+1);% 隐含层参数初始化
for i=1:num_net.num_hid2
    w1(i,:)=w1(i,:)/norm(w1(i,:));
end
w_o=rand(1,num_net.num_hid2);w_o=w_o/norm(w_o);% 输出层参数初始化
h=@(x,w)1/(1+exp(-sum(w.*x)));% sigmoid神经元函数
network_fullconnect=struct('num_net',num_net,'w_hid1', ...
    w1,'w_out',w_o,'Actfun',h);% 声明神经网络，其中Actfun为激活函数
% 利用梯度下降法迭代更新系数：
mn=1000;% 最大迭代次数
alpha=0.01;% 学习率
accu=zeros(mn,1);
for k=1:mn
    tic
    accu0=0;
    for i=1:n
        ho=zeros(1,network_fullconnect.num_net.num_hid2);
        for i1=1:network_fullconnect.num_net.num_hid2
            ho0=network_fullconnect.Actfun(x(i,:),network_fullconnect.w_hid1(i1,:));
            ho(i1)=ho0;
        end
        y_pre=network_fullconnect.Actfun(ho,network_fullconnect.w_out);
        w_o0=(y_pre-y(i))*ho*exp(-sum(network_fullconnect.w_out.*ho))/...
            (1+exp(-sum(network_fullconnect.w_out.*ho)))^2;
        network_fullconnect.w_out=network_fullconnect.w_out-alpha*w_o0;% 修正输出层系数
        for i2=1:network_fullconnect.num_net.num_hid2
            w10=(y_pre-y(i))*w_o0(i2)*x(i,:)*...
            exp(-sum(network_fullconnect.w_hid1(i2,:).*x(i,:)))/...
                (1+exp(-sum(network_fullconnect.w_hid1(i2,:).*x(i,:))))^2;
            network_fullconnect.w_hid1(i2,:)=network_fullconnect.w_hid1(i2,:)-alpha*w10;
        end
        accu0=accu0+abs(y_pre-y(i))/y(i);
    end
    accu(k)=accu0/mn;
    toc
end
% 可以看出基本的全连接神经网络具有一定运算精度。随着迭代次数的增加，精度不断
% 提高。
% 对于本节数据，本节介绍matlab神经网络预测设计的两种方案：
% (1)通过回归学习器和分类学习器可以进行全连接神经网络预测模型设计，对于本节
% 数据可见函数文件：trainRegressionModel_ANN.m，或用回归学习器打开：
% RegressionLearnerSession_ANN。
% (2)通过神经网络拟合可以设计具有sigmoid激活函数的隐含层和线性函数的输出层的
% 全连接神经网络。当本节数据在工作区时，可打开文件：Neural_network_fitting.m
% 进行训练验证数据，预测数据用y=net(x)进行。
% 利用matlab神经网络函数设计全连接神经网络模型：该部分可见文件：
% Fully_connected_neural_network.m。
input=x;% 输入数据，列数代表样本数
output=y;% 响应
n=length(output);
n1=ceil(0.7*n);% 训练数据个数
% 将指标变为列向量，划分训练样本和测试样本
input_train=input(1:n1,:)';
output_train=output(1:n1,:)';
input_test=input(n1+1:end,:)';
output_test=output(n1+1:end,:)';
% 节点数设置
inputnum=10;% 输入层节点数，即x维度数
hiddennum=10;% 隐含层数目
outputnum=1;% 输出层数目
% 数据归一化处理
[inputn,inputps]=mapminmax(input_train);% 归一化到[-1,1]之间
[outputn,outputps]=mapminmax(output_train);
% 构建网络
net=newff(inputn,outputn,hiddennum,{'tansig','purelin'},'trainlm');% 模型
% 隐含层使用tansig激活函数，输出层使用purelin线性函数，使用梯度下降法训练
W1=net.IW{1,1};% 输入层到中间层权值
B1=net.b{1};% 中间各层神经元阈值
W2=net.LW{2,1};% 中间层到输出层权值
B2=net.b{2};% 输出层各神经元阈值
% 训练参数配置
net.trainParam.epochs=1000;% 训练迭代次数
net=train(net,inputn,outputn);% 进行训练
% 测试
inputn_test=mapminmax('apply',input_test,inputps);
an=sim(net,inputn_test);
test_simu=mapminmax('reverse',an,outputps);% 对预测结果反归一化
%% 2.卷积神经网络Convolutional neural network
% 卷积神经网络在全连接神经网络基础上增加了卷积计算和池化层，对输入数据进行了
% 特征提取。
% 卷积计算层将输入数据与滤波器进行卷积，得到了输入数据的特征量。卷积神经网络
% 一般应用在图像识别邻域，因此本节说明中输入为一矩阵。在进行卷积计算时，需要
% 声明：滤波器个数，步长，填充值。每个滤波器进行卷积后得到该滤波器下的输出。
% 计算完卷积得到输出后，利用ReLU激活函数计算激励，然后送入池化层，池化层对于
% 输入数据取某一区域内平均或最大。
% 本节使用matlab内置图片集数据学习训练神经网络
clear,clc
dataFolder = fullfile('DigitDataset');% 将图片存入文件夹中，一类图片放在一个
% 子文件夹内
imds = imageDatastore(dataFolder, ...
    'IncludeSubfolders',true, ...
    'LabelSource','foldernames');% 读取图片，以图片所在文件夹名字作为图片标签
% CNN可以由深度网络设计器进行设计，可见文件mlx.CNN.mlx，将文件与mat文件：
% params_2022_12_25__21_33_12.mat放在同一文件夹中进行调用，其主要建模步骤如
% 下所示：
% 加载参数项：mat文件
trainingSetup = load("params_2022_12_25__21_33_12.mat");
% 导入数据
% 导入训练和验证数据。
imdsTrain = trainingSetup.imdsTrain;
[imdsTrain, imdsValidation] = splitEachLabel(imdsTrain,0.7);
% 调整图像大小以匹配网络输入层。
augimdsTrain = augmentedImageDatastore([227 227 1],imdsTrain);
augimdsValidation = augmentedImageDatastore([227 227 1],imdsValidation);
% 设置训练选项
% 指定训练时要使用的选项。
opts = trainingOptions("sgdm",...
    "ExecutionEnvironment","auto",...
    "InitialLearnRate",0.01,...
    "Shuffle","every-epoch",...
    "Plots","training-progress",...
    "ValidationData",augimdsValidation);
% 创建层组
layers = [
    imageInputLayer([227 227 1],"Name","imageinput","Normalization","zscore")
    convolution2dLayer([5 5],6,"Name","conv","Padding","same")
    reluLayer("Name","relu")
    maxPooling2dLayer([5 5],"Name","maxpool","Padding","same","Stride",[2 2])
    convolution2dLayer([5 5],16,"Name","conv_1","Padding","same")
    reluLayer("Name","relu_1")
    maxPooling2dLayer([5 5],"Name","maxpool_1","Padding","same","Stride",[2 2])
    convolution2dLayer([5 5],12,"Name","conv_2","Padding","same")
    reluLayer("Name","relu_2")
    maxPooling2dLayer([5 5],"Name","maxpool_2","Padding","same","Stride",[2 2])
    fullyConnectedLayer(30,"Name","fc")
    reluLayer("Name","relu_3")
    fullyConnectedLayer(10,"Name","fc_1")
    softmaxLayer("Name","softmax")
    classificationLayer("Name","classoutput")];
% 训练网络
% 使用指定的选项和训练数据对网络进行训练。
[net, traininfo] = trainNetwork(augimdsTrain,layers,opts);
%% 3.循环神经网络recurrent neural network
% 卷积神经网络卷积层从输入空间维度提取特征，而循环神经网络从时间维度上处理数
% 据。相比较全连接神经网络，循环神经网络与卷积神经网络卷积层相同，首先进行权
% 值共享，减少参数数量，然后增加一为0的存储单元h和其对应的权值，按照输入序列
% 的顺序，进行循环。首先与输入序列中第一项做运算，然后把该运算结果和神经元运
% 算结果相加作为输入送入第二神经元内，和第二项输入一同进行运算，以此类推。
% RNN循环过程可表示为：s(t)=f(w1*x(t)+w2*s(t-1));o(t)=g(w1*s(t))，其中s(t)为
% 循环输出项，o(t)为神经元输出。
% 本节通过GRU和LSTM两种RNN神经网络基础设计方案进行数据预测。本节训练数据选自
% 文件负荷天气数据中负荷前10000项，通过函数文件data_process.m进行数据输入输
% 出划分，具体程序如下：
clc;clear;close all
% 读取数据预处理
data=readtable('负荷天气数据.csv','range','B2:B10001');
data=table2array(data);
[x,y]=data_process(data,30);%前30个时刻 预测下一个时刻
% 归一化
[xs,mappingx]=mapminmax(x',0,1);x=xs';
[ys,mappingy]=mapminmax(y',0,1);y=ys';
% 划分数据
n=size(x,1);
m=round(n*0.7);%前70%训练，对最后30%进行预测
XTrain=x(1:m,:)';
XTest=x(m+1:end,:)';
YTrain=y(1:m,:)';
YTest=y(m+1:end,:)';
% 设计网络
numFeatures = size(XTrain,1);% 输入特征数
numResponses = 1;% 输出响应数
numHiddenUnits = 200;% 隐含层数目
%  GRU设计
layers1 = [
    sequenceInputLayer(numFeatures,"Name","sequence")
    gruLayer(numHiddenUnits,"Name","gru")
    fullyConnectedLayer(numResponses,"Name","fc")
    regressionLayer("Name","regressionoutput")];
% LSTM设计
layers2 = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    fullyConnectedLayer(numResponses)
    regressionLayer];
% 训练优化设计
options = trainingOptions('adam', ...% adam形式梯度下降优化
    'MaxEpochs',200, ...% 最大迭代次数
    'GradientThreshold',1, ...
    'InitialLearnRate',0.001, ...% 初始学习率
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',125, ...
    'LearnRateDropFactor',0.2, ...
    'Verbose',0);
% 训练两种方案网络
net1 = trainNetwork(XTrain,YTrain,layers1,options);
net2 = trainNetwork(XTrain,YTrain,layers2,options);
numTimeStepsTest = size(XTest,2);
% 对网络1结果分析
for i = 1:numTimeStepsTest
    [net1,YPred(:,i)] = predictAndUpdateState(net1,XTest(:,i),'ExecutionEnvironment','cpu');
end
% 反归一化
predict_value=mapminmax('reverse',YPred,mappingy);
true_value=mapminmax('reverse',YTest,mappingy);
disp('结果分析')
rmse=sqrt(mean((true_value-predict_value).^2));
disp(['根均方差(RMSE)：',num2str(rmse)])
mae=mean(abs(true_value-predict_value));
disp(['平均绝对误差（MAE）：',num2str(mae)])
mape=mean(abs(true_value-predict_value)/true_value);
disp(['平均相对百分误差（MAPE）：',num2str(mape*100),'%'])
fprintf('\n')
figure
plot(true_value,'-*','linewidth',1)
hold on
plot(predict_value,'-s','linewidth',1)
legend('实际值','预测值')
grid on
title('GRU')
% matlab深度网络设计：
% 通过matlab深度网络设计app可以设计不同形式的神经网络进行训练预测，下面简单
% 对设计器中各层作用进行介绍。
% matlab神经网络设计使用一般按照如下形式进行：
% (1). layers=[....];% 设计网络中所包含的层
% (2). options = trainingOptions(....);% 设计求解训练的优化策略及参数
% (3). net = trainNetwork(XTrain,YTrain,layers,options);% 导入训练数据进行训练
% (4).[net1,YPred(:,i)] = predictAndUpdateState(net1,XTest(:,i),...
% 'ExecutionEnvironment','cpu');% 进行预测并更新模型，声明执行环境
% 本节介绍layers部分设计：
tempLayers1 = [
    sequenceInputLayer(8,"Name","sequence","Normalization","zscore")% 序列
    % 形式输入，输入个数为8，对输入进行归一化处理
    convolution1dLayer(3,32,"Name","conv1d","Padding","same")% 一维卷积计算
    % 层，卷积核长度为3，有32个卷积核
    reluLayer("Name","relu")% 激活函数为relu
    maxPooling1dLayer(5,"Name","maxpool1d","Padding","same")% 最大池化层，每
    % 5项取最大
    lstmLayer(128,"Name","lstm")% lstm网络层，有128个单元
    fullyConnectedLayer(10,"Name","fc")% 全连接层
    regressionLayer("Name","regressionoutput")];% 回归输出
tempLayers2 = [
    featureInputLayer(1,"Name","featureinput","Normalization","zscore")% 特
    % 征形式输入
    fullyConnectedLayer(10,"Name","fc_1")% 全连接层
    tanhLayer("Name","tanh")% 激活函数为tanh
    classificationLayer("Name","classoutput_1")];% 分类输出
tempLayers3 = [
    imageInputLayer([227 227 3],"Name","imageinput")% 输入图像，大小为：
    % 227*227*3
    convolution2dLayer([3 3],32,"Name","conv","Padding","same")% 二维卷积层
    reluLayer("Name","relu_1")
    maxPooling2dLayer([5 5],"Name","maxpool","Padding","same")
    fullyConnectedLayer(10,"Name","fc_2")
    classificationLayer("Name","classoutput")];
%% 4.智能元胞自动机Intelligent cellular automata
% 本节通过设计生命游戏构建一个智能元胞自动机模型。
% 现有一座小岛，岛上部分资源有限，而部分资源随机产生，岛上有森林，河流，老虎
% ，兔子，老鹰。现有一批轮船失事幸存者流落在该岛上。你需要通过各资源及条件限
% 制，模拟小岛的动态变化过程。
% 有如下任务需要幸存者完成：
% (1)建造避难所：用来过夜，逃避老虎袭击。每次幸存者与老虎相遇，若不能及时逃回
% 避难所，将被老虎吃掉1人；每次过夜，未在避难所者将减少10%，至少1人。建造避难
% 所需要人发挥力量属性完成。
% 避难所建造条件：当避难所数量不够人住时，在总水源获取达50%和能量50%时，首要
% 任务为建造避难所。建造避难所需要消耗一定的能量、水源和树木，且建造地点只能
% 在森林。
% (2)获取水源：用来补充幸存者水资源，需要力量属性，消耗一定的能量。获取水源只
% 能在河流获取。
% (3)获取能量：用来补充能量，需要力量属性和智慧属性协同作用，需要消耗一定的水
% 源。获取能量的方式为捕食兔子。
% (4)附加限制：森林会不断更新，老虎会捕食兔子，老鹰也会捕食兔子，同样老虎也会
% 杀害幸存者。当某避难所内既有男性又有女性时，有几率发生繁衍行为。
% (@)幸存者一共有10人，5人为女性，5人为男性，每个人都具有力量，智慧两项属性，
% 在执行任务时，每个人只能发挥一项属性。力量能力与利用森林资源建造避难所，获
% 取水源两项任务直接关联；智慧能力可以完成逃避老虎追击，捕食兔子的任务，但需
% 要力量能力协助完成。每个人都有能量，水源两项需求，当满足这两项需求到一定程
% 度时，才能开始各项任务，当任意一项需求为0持续1天时，该幸存者将会死亡。
clear,clc
% 幸存者初始化
num=10;
gender=[zeros(1,5),ones(1,5)];% 性别：0代表女性，1代表男性
locationx=cumsum(ones(1,10));% 坐标x
locationy=1e3*ones(1,10);% 坐标y
power=80+20*rand(1,10);% 能量
wisdom=80+20*rand(1,10);% 智慧
water=70+30*ones(1,10);% 水源
energy=70+30*ones(1,10);% 能量
sur=struct('id',locationx,'gender',gender,'location', ...
    [locationx;locationy],'power', ...
    power,'wisdom',wisdom,'water',water,'energy',energy);
% 兔子位置初始化
num1=5000;
rabbit=1e3*rand(2,num1);
% 老虎位置初始化
num2=5;
tiger=400+300*rand(2,num2);
% 老鹰位置初始化
num3=100;
eagle=1e3*rand(2,num3);
% 场地初始化
x=1e3;y=1e3;
field=zeros(x,y);
for i1=1:x
    for i2=1:y
        if(i2<=i1+300&&i2>=i1-400)
            field(i1,i2)=5;% 令3，4，5代表树木，数字不同代表茂密程度不同
        end
        if(i2<=300&&i2>=230&&i1<=535)
            field(i1,i2)=1;% 令1代表河流
        end
        if(i1<=535&&i1>=465&&i2<=770&&i2>=230)
            field(i1,i2)=1;% 令1代表河流
        end
        if(i2<=770&&i2>=700&&i1>=465)
            field(i1,i2)=1;% 令1代表河流
        end
    end
end
field=repmat(field,2,2);
% 避难所初始化
shelter=zeros(3,1);lim=4;% 一个避难所容量为4
% 迭代计算
mn=50;% 迭代总时长
for i=1:mn
    tic
    if(mean(sur.energy)>=50&&mean(sur.water)>=50)
        if(lim*sum(shelter(1,:))<num)
            meanloc=[mean(sur.location(1,:));mean(sur.location(2,:))];
            tod=300+400*rand(2,1);tod=(tod-meanloc)/norm(tod-meanloc);
            alpha=sum(sur.energy)/(num*100);alpha=400*alpha;
            loc2=floor(meanloc+alpha*tod);
            for k0=1:2
                if(loc2(k0)>=x)
                    loc2(k0)=x-1;
                elseif(loc2(k0)<=0)
                    loc2(k0)=2;
                end
            end
            sur.energy=sur.energy-10;sur.water=sur.water-10;
            if(field(loc2(1),loc2(2))>=3&&field(loc2(1),loc2(2))<=5)
                shelter=[shelter,[1;loc2(1);loc2(2)]];
                field(loc2(1),loc2(2))=10;% 10代表避难所
                field(loc2(1)-1,loc2(2))=0;% 建造一所避难所产生一块空地
            end
            sur.location=repmat(loc2,1,num);
        else
            sur.energy=sur.energy-10;sur.water=sur.water-10;
            num=num+1;% 当能量，水源，避难所数量都足够时，增加人口数量
            sur.gender=[sur.gender,floor(rand())];
            sur.id=[sur.id,sur.id(end)+1];
            locationx0=shelter(2,end);
            locationy0=shelter(3,end);
            sur.location=[sur.location,[locationx0;locationy0]];
            power0=50+20*rand();
            sur.power=[sur.power,power0];
            wisdom0=50+20*rand();
            sur.wisdom=[sur.wisdom,wisdom0];
            water0=50+30*rand();
            sur.water=[sur.water,water0];
            energy0=50+30*rand();
            sur.energy=[sur.energy,energy0];
        end
    elseif(mean(sur.energy)<50)
        for k1=1:num
            sur.location(:,k1)=sur.location(:,k1)+-300+300*rand(2,1);
            sur.energy(k1)=sur.energy(k1)-5;
            sur.water(k1)=sur.water(k1)-5;
        end
        for k1=1:num
            po=sur.location(:,k1);
            d=[];dn=0;dr=0;dt=0;de=0;
            for k2=1:num % 每个人身边队友数量
                if(k2~=k1&&sqrt(sum((sur.location(:,k2)-po).^2))<=100)
                    d=[d,[k2;sqrt(sum((sur.location(:,k2)-po).^2))]];
                    dn=dn+1;
                end
            end
            for k2=1:num1 % 遇到兔子的数量
                if(sqrt(sum((rabbit(:,k2)-po).^2))<=30)
                    dr=dr+1;
                end
            end
            for k2=1:num2 % 遇到老虎的数量
                if(sqrt(sum((tiger(:,k2)-po).^2))<=30)
                    dt=dt+1;
                end
            end
            for k2=1:num3 % 遇到老鹰的数量
                if(sqrt(sum((eagle(:,k2)-po).^2))<=70)
                    de=de+1;
                end
            end
            tm=size(d);
            if(tm(1)~=0)
                if(dt>=1)% 遇到老虎时处理
                    wis0=0;pow0=0;
                    for k3=1:length(d(1,:))
                        wis0=wis0+sur.wisdom(d(1,k3));
                        pow0=pow0+sur.power(d(1,k3));
                    end
                    if((2/3*wis0+1/3*pow0)/(dn+1)>90)
                        sur.location(:,k1)=sur.location(:,k1)+-200+200*rand(2,1);
                        for k3=1:length(d(1,:))
                            sur.location(:,d(1,k3))=sur.location(:,d(1,k3))+...
                                -200+200*rand(2,1);
                        end
                    else
                        sur.location(:,k1)=[-1e3;-1e3];
                        for k3=1:length(d(1,:))
                            sur.location(:,d(1,k3))=sur.location(:,d(1,k3))+...
                                -200+200*rand(2,1);
                        end
                    end
                else
                    sur.energy(k1)=sur.energy(k1)+20*(dr-de)/(dn+1);
                end
            else
                sur.location(:,k1)=[-1e3;-1e3];
            end
        end
        for k1=num:-1:1
            if(sur.energy(k1)<=0||sur.location(1,k1)<=-900||sur.water(k1)<=0)
                num=num-1;
                sur.energy(k1)=[];
                sur.gender(k1)=[];
                sur.id(k1)=[];
                sur.location(:,k1)=[];
                sur.power(k1)=[];
                sur.water(k1)=[];
                sur.wisdom(k1)=[];
            end
        end
        for k1=1:num
            if(sur.energy(k1)>=100)
                sur.energy(k1)=100;
            end
            for k2=1:2
                if(sur.location(k2,k1)>=x)
                    sur.location(k2,k1)=x-1;
                elseif(sur.location(k2,k1)<=1)
                    sur.location(k2,k1)=2;
                end
            end
            if(sur.water(k1)>=100)
                sur.water(k1)=100;
            end
        end
    elseif(mean(sur.water)<50)
        meanloc=[mean(sur.location(1,:));mean(sur.location(2,:))];
        tod=300+400*rand(2,1);tod=(tod-meanloc)/norm(tod-meanloc);
        alpha=sum(sur.energy)/(num*100);alpha=400*alpha;
         for k1=1:num
             sur.location(:,k1)=sur.location(:,k1)+alpha*tod;
             sur.energy(k1)=sur.energy(k1)-5;
             sur.water(k1)=sur.water(k1)-5;
         end
         for k1=1:num
             for k2=1:2
                 if(sur.location(k2,k1)>=x)
                     sur.location(k2,k1)=x-1;
                 elseif(sur.location(k2,k1)<=0)
                     sur.location(k2,k1)=2;
                 end
             end
         end
        for k1=1:num
            po=sur.location(:,k1);
            d=[];dn=0;dt=0;
            for k2=1:num % 每个人身边队友数量
                if(k2~=k1&&sqrt(sum((sur.location(:,k2)-po).^2))<=100)
                    d=[d,[k2;sqrt(sum((sur.location(:,k2)-po).^2))]];
                    dn=dn+1;
                end
            end
            for k2=1:num2 % 遇到老虎的数量
                if(sqrt(sum((tiger(:,k2)-po).^2))<=30)
                    dt=dt+1;
                end
            end
            tm=size(d);
            if(tm(1)~=0)
                if(dt>=1)% 遇到老虎时处理
                    wis0=0;pow0=0;
                    for k3=1:length(d(1,:))
                        wis0=wis0+sur.wisdom(d(1,k3));
                        pow0=pow0+sur.power(d(1,k3));
                    end
                    if((2/3*wis0+1/3*pow0)/(dn+1)>90)
                        sur.location(:,k1)=sur.location(:,k1)+-200+200*rand(2,1);
                        for k3=1:length(d(1,:))
                            sur.location(:,d(1,k3))=sur.location(:,d(1,k3))+...
                                -200+200*rand(2,1);
                        end
                    else
                        sur.location(:,k1)=[-1e3;-1e3];
                        for k3=1:length(d(1,:))
                            sur.location(:,d(1,k3))=sur.location(:,d(1,k3))+...
                                -200+200*rand(2,1);
                        end
                    end
                elseif(field(floor(sur.location(1,k1)),floor(sur.location(2,k1)))==1)
                    pow0=0;
                    for k3=1:length(d(1,:))
                        pow0=pow0+sur.power(d(1,k3));
                    end
                    for k3=1:length(d(1,:))
                        sur.water(d(1,k3))=sur.water(d(1,k3))+20*pow0/(100*(dn+1));
                    end
                end
            else
                sur.location(:,k1)=[-1e3;-1e3];
            end
        end
        k1=num;
        while k1>=1
            if(sur.energy(k1)<=0||sur.location(1,k1)<=-900||sur.water(k1)<=0)
                num=num-1;
                sur.energy(k1)=[];
                sur.gender(k1)=[];
                sur.id(k1)=[];
                sur.location(:,k1)=[];
                sur.power(k1)=[];
                sur.water(k1)=[];
                sur.wisdom(k1)=[];
            end
            k1=k1-1;
        end
        for k1=1:num
            if(sur.energy(k1)>=100)
                sur.energy(k1)=100;
            end
            for k2=1:2
                if(sur.location(k2,k1)>=x)
                    sur.location(k2,k1)=x-1;
                elseif(sur.location(k2,k1)<=1)
                    sur.location(k2,k1)=2;
                end
            end
            if(sur.water(k1)>=100)
                sur.water(k1)=100;
            end
        end
    end
    toc
    disp(num2str(num));
    pause(0.3)
    subplot(2,2,1)
    scatter(sur.location(1,:),sur.location(2,:))
    subplot(2,2,2)
    bar(sur.energy)
    subplot(2,2,3)
    bar(sur.water)
    km=rand();
    if(sum(shelter(1,:))~=0)
        if(km<=0.05)
            shelter(:,end)=[];% 避难所有几率被销毁
        end
    end
    % 兔子位置更新
    rabbit=1e3*rand(2,num1);
    % 老虎位置更新
    tiger=1e3*rand(2,num2);
    % 老鹰位置更新
    eagle=1e3*rand(2,num3);
    if(num==0)
        break
    end
end
% 通过改变field场地设置或令field动态变化，即可进行元胞自动机生命游戏设计。
